An efficient parallel method for relaxing to the minimum action wavefunction 
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Efficient and accurate numerical propagation of the time dependent Schrodinger equation is a 
problem with applications across a wide range of physics. This paper develops an efficient, trivially 
parallelizeable method for relaxing a trial wavefunction toward a variationally optimum propagated 
wavefunction which minimizes the propagation error relative to a platonic wavefunction which obeys 
the time dependent Schrodinger equation exactly. This method is shown to be well suited for 
incorporation with multigrid methods, yielding rapid convergence to a minimum action solution 
even for Hamiltonians which are not positive definite. 



As the general wave equation for nonrelativistic quan- 
tum mechanics, the time dependent Schrodinger equation 
holds sway over a broad range of modern physics. How- 
ever, due in part to its generality, the TDSE remains a 
challenging problem to treat numerically. For a wave- 
function expanded in a basis set, with operators repre- 
sented by matrices, the wide range of physical problems 
described by the TDSE makes it difficult to define prop- 
agation procedures which are both general and efficient. 
Nearly any simplifying assumption about the form of the 
Hamiltonian matrix will be violated by some problem of 
interest, while procedures which do not make such as- 
sumptions may be computationally expensive. 

The need for an efficient, general propagator for the 
TDSE has been highlighted in recent years by the rise of 
strong field physics, where problems often have partic- 
ularly poor numerical properties. A typical strong field 
experiment may involve an electron tunneling free from 
a molecule due to the field of an intense laser. After 
reaching the continuum, the electron is accelerated by 
the field of the laser and may traverse hundreds of bohr 
before returning to scatter energetically from the parent 
ion. The electric field due to the laser and the contin- 
uum electron may excite the ion, while a sufficiently en- 
ergetic recollision may liberate additional electrons. A 
strong field problem may thus include singular and irreg- 
ularly shaped molecular potentials, large length scales, 
high energies, and time dependent, nonperturbative ex- 
ternal forces. At present, calculations involving a single 
electron in three dimensions refiects a significant com- 
putational challenge [l|, while two electron calculations 
may require thousands of processors [2!] . 

Recently, the problem of propagating a time dependent 
wavefunction has been approached from the perspective 
of minimizing the action accumulated when the wave- 
function is expanded in a particular basis Q . For a finite 
basis set, this action is not necessarily zero. Minimizing 
this action was shown to minimize the deviation between 
the propagated wavefunction and a platonic "true" wave- 
function which obeys the time dependent Schrodinger 
equation exactly. By appropriate choice of basis set, the 
propagated wavefunction can be made to converge expo- 
nentially to the true wavefunction. The least action prop- 



agator thus represents a general method for accurately 
propagating the time dependent Schrodinger equation. 
However, the accuracy of this propagator is purchased at 
the cost of solving a potentially very large linear system 
of equations at every time step. For problems involving 
large or unstructured basis sets, where efficient propaga- 
tors are most needed, direct solution of this linear system 
is likely to prove prohibitively expensive. 

This paper addresses the problem of solving the least 
action linear system by developing an iterative method 
for decreasing the action accumulated by a trial wave- 
function. By dividing a spacetime volume into many 
small subvolumes, the accumulated action is decreased 
by finding corrections which minimize the action accu- 
mulated in each subvolume. For the common case of a 
local Hamiltonian, corrections to the trial wavefunction 
in nonoverlapping subvolumes are independent of one 
another, making this relaxation procedure trivially par- 
allelizeable. Combining this relaxation procedure with 
multigrid algorithms yields an overall solution procedure 
which replaces the solution of a single large linear system 
with the solution of many small linear systems, greatly 
reducing the computational cost. As multigrid methods 
have traditionally been restricted to use with linear sys- 
tems arising from elliptical PDEs, this reflects an expan- 
sion of the class of problems for which these powerful 
methods are suitable. 



I. REVIEW OF THE LEAST ACTION 
PRINCIPLE 

The least action propagator was derived in Q as a 
way to minimize the deviation between the propagated 
wavefunction and a platonic "true" wavefunction which 
obeys the time dependent Schrodinger equation exactly. 
While [3] employed a Taylor series expansion for this pur- 
pose, a simpler derivation requires only the fundamental 
theorem of calculus. If the true wavefunction is given 
by 'ip{x, t) and the approximate wavefunction is given by 
0(a;, i) — ^ CinXi{^)Tn{t) , the approximation error is 
given by S{x, t) = ijj{x, t) — (t>{x, t). The norm of this de- 
viation at some time t for a propagation step beginning 
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at ^0 is given by 



{s{t)m)={sitomto)) 



dt'—{5{t'mt')). 



The assumption that (i^ — H) \^) — can be used to 
ehminate all terms involving or \^) from the second 
term, so that 



{5{t)\5{t)) = {5{tomh)) - 2z / dt' {5[t')\ H \5{t')) 

(2) 

+2ij dt'm')\i^j^-H)\m)- 



S + E^^^f*, where fi = C,nTn{t) - {x^\i^{x,t)), min- 
imizing 5" with the constraint that fi = for all i yields 
the minimum action solution with the initial condition 

When using the calculus of variations in this way, it 
must be remembered that it is, strictly speaking, impos- 
sible to "minimize" a complex quantity such as the ac- 
tion. This difficulty is particularly relevant in the context 
of seeking an iterative solution procedure, as the conver- 
gence of such a procedure can be more easily judged in 
the context of minimizing a real quantity than in ap- 
proaching an extremum of a functional derivative. Ac- 
cordingly, it is useful to define an action residual 



In this equation, terms involving (j>(x, t) result from 
imperfectly propagating the wavefunction, while terms 
involving 5{x, t) result from the inability to represent the 
true wavefunction in the chosen basis. As the propa- 
gation error is proportional to the integral of the La- 
grangian density L = i-^ — H, the wavefunction which 
minimizes propagation error is that which minimizes the 
action accumulated over the time step. 

The least action condition can be turned into a prop- 
agation scheme using the calculus of variations. For a 
Hamiltonian H — Hq + V{x,t) with matrix elements in 
space and time given by 



H, 



t+At 



dxx*ix)HoXjix) 



(3) 



V.jnm = I dt' dxX*{x)T*{t)V{x,t)Xj{x)Tra{t) 



Ot.j = / dxx*{x)xA^) 



(4) 
(5) 
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as a quantity to measure the deviation of a wavefunction 
from the action-extremizing condition given by Eq. 1^1 
Using this definition, the norm of the residual 
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is positive definite, and 
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(11) 



(12) 



when the residual is zero and Eq. IHlis satisfied. Restating 
the least action problem in this way also gives a simple 
path toward improving a trial wavefunction. If a trial 
wavefunction gives a nonzero residual in Eq. ^ the ex- 
pansion coefficients must be adjusted so as to decrease 
the overall norm of the action residual. Despite their 
technical inaccuracy, the terms "least action" and "min- 
imizing the action" will be retained, as the terminology 
is familiar from classical physics. 



dtT:{t)Trn{t) 



/ dtT* {t)T„,{t), 



(6) 



(7) 



the action accumulated between t and t + At is given by 
<l}{x,t) = J2in'^inXi{x)Tn{t) IS givcu by 

i.j,n.m 

and is minimized when 
dS 



for all j,m. Initial and boundary conditions can be 
specified using Lagrange multipliers. Defining S' = 



II. ITERATIVE RELAXATION TO THE LEAST 
ACTION SOLUTION 



The least action propagator purchases its rapid con- 
vergence to the true wavefunction at a high computa- 
tional cost. For a basis with A^^; spatial and Nt temporal 
basis functions, the least action propagator must solve 
for Nx{Nt + 1) coefficients (including Lagrange multipli- 
ers) at every time step. For the large and irregular spa- 
tial bases where efficient propagation schemes are most 
needed, the cost of directly solving the least action linear 
system defined by Eq. IHlis likely to prove unacceptably 
high. 

As an alternative to direct solution, many iterative al- 
gorithms have been developed for the solution of linear 
systems of equations. Methods such as the conjugate 
gradient or generalized minimum residual algorithm at- 
tempt to solve a general linear system Ax = & by solving 
a smaller linear system within a Krylov subspaceQ. 



3 



A common theme in such iterative methods is the need 
for a preconditioner M^^ which has the effect of approx- 
imately inverting the matrix A. Given a trial solution x 
with residual r = b — Ax, a correction term 5x = M~^r 
can be used to improve the trial solution, so that a new 
trial solution is given hy x' — x -\- 6x. In the limit that 
= A~^, the correction yields an exact solution; how- 
ever, it is more typical to choose a preconditioner which 
is simple to calculate and yields acceptable performance. 
Generic preconditioners include Jacobi or Gauss Seidel 
relaxation, while physical knowledge may allow construc- 
tion of preconditioners well suited to individual problems. 

To define a preconditioner for the least action propa- 
gator, it is necessary to return to the action integral over 
a volume V 



dt / dV(f)*L<p, 



(13) 



where L — i-^ — H \s the Lagrangian operator and the 
Hamiltonian H = T + V is the sum of a kinetic energy 
operator T = and a (possibly time dependent) 

potential energy operator V . The divergence theorem 
can then be used to express the kinetic energy operator 
in terms of singly diflerentiable quantities 

— [ dt f dV(t)*V^(t)=— I dt I (j)*V^-dA+ 
2to J Jy '2m J Jqy 

dt I dNV(l>* ■ Vcj) ■ dA, 



1 

2m 



(14) 



allowing consideration of trial wavefunctions with discon- 
tinuous derivatives. Given a Lagrangian which is a local 
function of the spatial coordinates and some subvolume 
U of the original volume, the action integral can be bro- 
ken up into the action accumulated inside U and that 
accumulated outside: 



where 



-Tr- 
im, 

1 



S = Sv + Sfj, 



(15) 



„ , dt / dUV(/)* • V0 • dA+ 
2m ./ Jo 



dt / di](j)*V{x,t)(l) 

u 



(16) 



and CF is the subvolume of V not contained in U. 

It can now be seen that the requirement that the prop- 
agated function (j) minimize the action accumulated over 
all space places restrictions on ^u, the function restricted 
to U. In order for (j) to minimize the action over all space, 
must minimize the action accumulated over U, sub- 
ject to the requirement that ^u|au = ^ul^u- If this were 
not the case, it would be possible to construct a relaxed 
wavefunction which would come closer to minimizing the 
action. 



Here it is convenient to consider the norm of the least 
action residual. Being positive definite, this norm allows 
two candidate wavefunctions to be compared by the de- 
gree to which they minimize the action. Let iV((/)) be 
the norm of the total action residual, with iVu((^u) and 
■^\}{4'\j) the norms of the action residuals inside and out- 
side the volume of interest. If there exists some function 
(fv for which Lpv\dv = <t>ts\dv and iVu(<Pu) < -^%(0u), then 
it is possible to construct a new global wavefunction (j)' , 
where 4>'^ = (pu and 4>'^ = (p^, such that 

N{4>') = N^{^) + N^{4>^) < N{^) = 7Vu(<Au) +7Vo(</,p). 

(17) 

This logic leads naturally to a relaxation procedure. 
Given a trial solution (p and a collection of (possibly 
overlapping) subvolumes U, the trial wavefunction may 
be systematically improved by minimizing the action 
accumulated in each subvolume, subject to the restric- 
tion that each relaxation step leave the wavefunction un- 
changed on the boundary and the exterior of the relax- 
ation volume. For a local Lagrangian, relaxation over 
two nonover lapping volumes is independent, allowing for 
parallel execution. 

To cast this procedure in the form of a preconditioner, 
a relaxation step must find a correction Scpv such that 
LuSipij = —LTj(t>, with initial conditions 6(f>{x, to) = and 
boundary conditions (5(/)|aiu = 0. 

Here it must be noted that the use of Lagrange mul- 
tipliers means that the relaxed solution cj) + S(j) is not 
guaranteed to be an extremum of the action, but rather 
a critical point. For basis functions which are nonzero 
at the initial time or on the boundaries of the relaxation 
volume, g^'? 7^ 0, due to the nonzero Lagrange multi- 

pliers associated with imposing the boundary conditions. 
For this reason, the relaxation procedure described here 
is not guaranteed to decrease the action accumulated in 
the relaxation volume. As will be seen, increasing the 
size of the basis set within the relaxation volume helps 
to minimize this effect. 



III. MULTIGRID RELAXATION 

The relaxation procedure introduced in the previous 
section has the side effect of introducing a length scale 
into the problem which does not originate with the PDE 
being solved. Because the relaxation procedure cannot 
change the wavefunction outside of the relaxation vol- 
ume, the size of the relaxation volume affects the rate of 
convergence. If the space is divided up into a set of re- 
laxation volumes of size V, the relaxation procedure will 
quickly eliminate "rough" error terms which accumulate 
action over volumes smaller than V, while error terms 
which accumulate action over volumes larger than V will 
be eliminated more slowly. This problem can be rectified 
somewhat by relaxing over larger volumes, at the cost of 
increasing the size of the linear system which must be 
solved. 
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Multigrid approaches [1, @ seek to eliminate this de- 
pendence on the relaxation volume by defining a treelike 
hierarchical basis set. At the base level, a single element 
is defined over the entire volume of interest. Within this 
element, functions are defined in terms of some small ba- 
sis set - here, low order polynomials. To describe details 
of a function which cannot be described by this small 
basis set, the element is subdivided into n smaller ele- 
ments, each with their own associated basis. Thus, the 
volume is divided into one element at the zeroth level, 
n elements at the first level, at the second level, and 
so on. This repeated subdivision produces a treelike ba- 
sis, with the parents of a volume element e consisting of 
those elements which contain e as a subvolume and the 
children consisting of those elements which are a subvol- 
ume of e. The leaves of the tree consist of those elements 
which have not been subdivided ~ ie, the finest level of 
the mesh. 

In defining this treelike basis set, there is no a pri- 
ori need to subdivide all volumes equally, and it may be 
convenient to give a compressed representation of some 
function by using a few coarse elements to describe the 
wavefunction in some region where it contains little de- 
tail, but many fine elements in some region where it os- 
cillates more rapidly. In this way, the multigrid basis 
set is closely related to wavelet decomposition, and offers 
similar possibilities for reduction of memory and storage 
requirements. 

The removal of the relaxation volume from the rate of 
convergence of the multigrid relaxation procedure is ac- 
complished by transferring the problem between different 
levels of the tree, so that relaxation may occur over many 
different length scales. If P^~^ is a projection operator 
mapping a function defined on the A'^th level to a func- 
tion defined on the (coarser) N — 1st level, and is 
an interpolation operator transferring functions the other 
way, so that P^~^I^_i is the identity on the — 1st 
level, and I^_j^P^~^ has eigenvalues of 1 for functions 
which can be expanded in the A^ — 1st level basis and 
otherwise, a linear equation defined on one level can be 
projected onto another level. For a linear system 



on the coarse grid and interpolating to the fine grid, so 
that 



(18) 



defined on the fine level, the projected problem on the 
coarse level is given by 

{Pl:l-'A%Pr-'){Pl^-'x'') ^ P^-'b^, (19) 

where P^^^i^/^_,. 



Given a trial solution with residual — b^ 

4^ 



A^x^ , a correction Sx^ can be found either by solving 



(20) 



on the fine grid or 



5x 



N 



tN /piV-i 4JV pJ iV-iwi piV-i iV /rjr,-. 



{P^-'A^^Pl^-'){P^~Hx^) = P^-V^ (21) 



The advantages of finding a correction on the coarse 
grid are twofold. First, there are fewer basis functions 
defined on the coarse grid, so that the overall problem is 
smaller. Second, the elements defined on the coarse grid 
are larger than those defined on the fine grid by a factor of 
n, so that a relaxation scheme such as the one described 
in Section |TT] will eliminate error terms which accumu- 
late action in volume nV rather than V. By recursively 
transferring the problem between different levels of the 
hierarchy, error terms acquiring action over any length 
scale can be made to converge, ideally yielding an overall 
procedure which is rapidly convergent with rates of con- 
vergence which do not depend on the size of the finest 
element. 

To date, multigrid approaches have primarily been 
used for problems arising from elliptical PDEs, due pri- 
marily to the lack of a suitable relaxation method to 
use for corrections at each level of the hierarchy. Re- 
laxation methods such as Gauss-Seidel require a Hamil- 
tonian which is either positive or negative definite^. If 
the eigenvalue spectrum spans 0, some error terms will 
diverge rather than converging. For elliptical problems, 
convergence is guaranteed for Gauss Seidel relaxation 
by the positive definite Hamiltonian, while non ellipti- 
cal problems such as the Schrodinger equation require 
that the error terms be eliminated in some other way, 
such as direct solution of the linear system at some level 
or developing an improved relaxation procedure. Multi- 
grid approaches have been used to solve the Schrodinger 
equation in General reviews of preconditioning 

methods are given by [ll|, [3 ' while [H, |l3] discuss pre- 
conditioners for closely related Helmholtz equation. 



IV. IMPLEMENTATION WITH LEGENDRE 
POLYNOMIALS 

The relaxation procedure described in this paper at- 
tempts to reduce the computational cost of performing 
a least action step by restricting the volume over which 
the action is minimized. Implicit in this restriction is the 
assumption that reducing the volume in this way will re- 
duce the number of basis set coefficients to be solved for 
- ie, that the chosen spatial basis consists of localized 
basis functions which are nonzero only in a restricted 
range. This property is true of many popular basis sets 
such as finite elements or the hierarchical multigrid bases 
described in the previous section but not for, e.g., Gaus- 
sian basis sets of the kind commonly used in quantum 
chemistry. For the purposes of testing the least action 
relaxation procedure, a suitable basis is provided by a 
set of low order Legendre polynomials. The Legendre 
basis has the advantage that many of the matrices in the 
least action equation are both sparse and analytically 
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calculable, with well behaved matrix elements. In one 
dimension, the range [—1,1] of the Legcndrc functions 
can be interpreted as a single finite clement, or the inter- 
val between two grid points in the hierarchical basis. In 
one spatial dimension, the overlap matrices O and U are 
given by the orthogonality relation 



/. 



^ dyPn{y)Pm{y) = 2n"+l ' 
while the Q matrix is given by 



(23) 



(24) 



if n > m and mod(n — to, 2) = 1, otherwise. Matrix 
elements of the kinetic energy operator are given by 



+ j\dyP'i{y)P'M 



-i{i + i) + {f+j) 

2m 



(25) 



when i > j and mod(? — j, 2) = 0, otherwise. These 
formulas make use of the identities Pn(l) = 1, P'li^) = 
^^^^ and P„{-y) = (-l)"P„(y). For the purposes 
of this paper, Vq is assumed to be constant throughout 
the volume of interest, so that matrix elements of the 
potential energy are proportional to the overlap matrix. 

Intergrid transfer operators are constructed by subdi- 
viding the [—1,1] range of the Legendre polynomials into 
two smaller elements, one ranging from [—1,0], the other 
from [0,1], and defining Legendre bases in the two subele- 

ments. Because P„(^-^) is an nth order polynomial over 
the range of both subelements, each basis function of the 
large element can be expanded as the sum of basis func- 
tions defined in the small elements, yielding interpolation 
matrices 



ij 



-L 



■/Pi{ 



-)PM) 



for the small element on the left and 



dy'P^{- 



1, 



PAv') 



(26) 



(27) 



for the small element on the right, which convert from the 
coarse to the fine basis. Projection matrices are simply 
the transpose of the interpolation matrices, normalized 
so that (P^-^) -f- P^-^)) • (/(-^^ + /(^)) = I is the identity 
matrix in the parent basis. 



,(L/H) 



j{L/R)^{L/R) 



rWf^iL) j(L) 
^ik '^kj ^ji 



^ik '^kj ^ji 



(28) 



In defining these matrices, the integrals of the dimen- 
sionless parameter y range from —1 to 1. For an el- 
ement "box" of size AxAt, these integrals acquire di- 
mension, with Tij oc 2/ Ax, Oij oc Ax/2, Qnm oc 2 and 



Unm oc At/2, and the least action equation becomes 

Cin[iOijQnm^X-TijUnm^-OijUnmVAxAt] = (Nj,m. 

(29) 

Defining Pmax = this equation can be recast as 

CiniiOijQnni - 2KTijUnm - fOijUnm] = OVj, TO, (30) 

where k = is the action due to kinetic energy 

acquired by a particle with momentum Pmax in time At 
and 1/ = VAt is the action acquired due to potential 
energy. 

Parameterizing the least action equation in this way 
helps to clarify the role of intergrid transfer operators in 
speeding convergence. If the size of the box is doubled, as 
might be seen in the transfer from a fine to a coarse grid, 
K — > k/4, while z/ is unaffected. Transferring from a fine 
to a coarse grid is thus isomorphic to propagating for a 
time At ^ At/4 in a potential V — > AV. As the problem 
is transferred from fine to coarse grids, the kinetic energy 
term will disappear, while transferring from coarse to fine 
grids will make this term dominate. If faster convergence 
is needed on a grid where the potential term dominates, 
it may be desirable to accomplish a single timestep in two 
steps of size At/2, which will have the efii'ect of mapping 
v — >■ Z//2 and k — >■ k/2. 

Because the Legendre polynomials are defined only 
within the confines of a single element box, it is nec- 
essary when working in the Legendre basis to enforce 
functional continuity at element boundaries. If (j)^^^ = 
Em C''^h'i^\x)Tn{t) is defined in the region from x—Ax 

to X and (/.(■^) = J2inC:^"^^x:i"H-^)Tnif.) is defined in the 
region from a; to a; -|- Ax, functional continuity is main- 
tained by requiring that 



(31) 



Boundary and initial conditions are specified in a sim- 
ilar way. Requiring (p{x,t) to equal f{x,t) at the initial 
time in some element a yields 

E ^in (-1)" = E fin (-l)"Vi, a. (32) 

n n 

Matching a function on the left boundary of some element 
a' yields 

Y.c\:^\-ir = Y.fln\-^r^n, (33) 

i i 

while matching the function on the right boundary of 
element a" yields 

E^L""^ = E/i«"^^- 



(34) 



Ensuring functional continuity between element a' on the 
left and a" on the right yields 



E^'^ = E(-i)^^ 



(a") 



Vn. 



(35) 
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When projecting a function f{x,t) into this basis, a 
projection which enforces continuity can be found by 
solving the hnear system 

C'ln^+T.^"=f^^^y^,n (36) 

n 

cif^+E(-l)^^" = /»^^^'" (37) 

T.(^in -J2i-^y(^^n^ =0^"- (38) 

i i 

where /^^ ^ J dx J dtf{x,t)x*''"\x)T*{t) is the overlap 
of the basis and the expanded function over element a. 

The least action relaxation step can be performed in 
one of two ways. To solve for the relaxed function di- 
rectly, it is necessary to solve the least action equa- 
tions with boundary conditions given by the trial wave- 
function on the boundary of the relaxation volume. If 
4>ix,t) = J2i n ^inXi{^)Tn{t) is the trial wavefunction, 
the relaxed wavefunction (j)'{x,t) = n^'inXi{^)Tn{t) 
is found by solving the linear system 

c'i"\iOlf''Qn,n ~ 2KTt/u^m " ^^Of/[/„„] = OVj, m, (3 

(39) 

with initial condition 



boundary conditions 



a. 



(40) 



boundary conditions 

T.^^^(-^'-J2C^ni-mn, (41) 



and 



and internal continuity enforced by 

Ec:^^ = E(-imfW 



(42) 



(43) 



Alternatively, the least action linear system can be cast in 
the form of a correction S(f){x,t) = .^5CinXi(x)Tn{t) 
to a trial function (j}{x^t). If 

= -Cl:\tOlf^Qn,n-2KTtfUr.„^-iyOtfUnrn]yj,m,(3 



(44) 

is the residual of the least action equations, 5C\'^^ can be 
found by solving 



iition 

E^4?^("i)" = o^*'"' 



with initial condition 



(45) 
(46) 



and 



E'^C'i^^^OVr., 



and internal continuity enforced by 

E^^»^=E(-ir^^»^Vn. 



(47) 



(48) 



(49) 



Restating the problem in the form of solving for a cor- 
rection to a trial wavefunction has the additional advan- 
tage that all terms in the least action linear system can 
be acted upon by intergrid transfer operators, making 
the problem susceptible to multigrid approaches. 



V. CONVERGENCE 

The relaxation procedure described in this paper seeks 
to solve for the minimum action wavefunction for a given 
initial condition by iteratively decreasing the action accu- 
mulated by a series of trial wavefunctions. Convergence 
to the true minimum action wavefunction is achieved 
when the action can no longer be lowered, and the resid- 
ual of the least action equation is zero. 

The convergence of the relaxed wavefunction to the 
minimum action solution was tested by expanding a trial 
wavefunction of the form f{x,t) = e*('^'^+"*) into a Leg- 
endre basis consisting of either one element ranging from 
[—1,1] or two elements ranging from [—1,0] and [0,1]. 
Convergence was measured by taking the ratio of the 
magnitude of the accumulated action before and after re- 
laxation, or alternatively by taking the ratio of the norm 
of the action residuals. Rates of convergence using one 
element are shown in Figures [T] and [21 while rates of con- 
vergence using two elements are shown in Figures |3] and 2] 
Because k and v may vary widely from problem to prob- 
lem, or within the same problem due to intergrid transfer, 
calculations were made for a wide range of both param- 
eters, including k >> 1, corresponding to large amounts 
of action due to high kinetic energy, and v < where the 
energy eigenvalue spectrum spanned zero. For all values 
of K and u tested, the relaxation procedure yielded rapid 
rates of convergence in the limit of sufficiently large basis 
sets. Rates of convergence were somewhat higher for cal- 
culations made with larger basis sets, although this must 
be balanced against the higher cost of an individual re- 
laxation step. More important was the requirement that 
the basis set yield sufficient free parameters to allow for 
convergence while satisfying + Nt{2 + — 1) con- 
straint equations, where iVg is the number of elements. 
For Nx = Nt = 3 the 9 free parameters are matched 
by 9 constraint equations for a single element calculation 
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while a two element calculation yields 18 free parame- 
ters and 12 constraints. The large number of constraint 
equations, with corresponding action costs in the form 
of Lagrange multipliers, means that rates of convergence 
may be slow or even diverging for small basis sets. Higher 
order basis sets yielded rapid convergence over nearly all 
of parameter space, with convergence slowing only in the 
vicinity of the —w = + v curve where the trial so- 
lution is already a good approximation to the minimum 
action wavefunction. 

Comparing these figures, it is apparent that relaxing 
over two elements simultaneously is very inefficient. As 
the cost of solving a linear system for n unknowns scales 
as n^, relaxing over 2 volumes simultaneously is approxi- 
mately 4 times as expensive as relaxing over each volume 
separately. However, comparing figures [1] and [3J or [5] and 
m shows that this added expense does not result in ap- 
preciably faster rates of convergence. Despite this ineffi- 
ciency, the inability of the relaxation procedure to change 
the wavefunction on the boundary of the relaxation vol- 
ume means that a convergent procedure must at some 
point relax over volumes larger than a single element. 

A more efficient way of relaxing over large volumes 
makes use of intergrid transfer operators. In Figures [S] 
and [HI relaxation over two adjacent volumes was accom- 
plished by first transferring the residual from two fine 
elements to one coarse element using Eq. The cor- 
rection found on the coarse grid was then transferred 
back to the fine grid using Eq. [211 This partially relaxed 
wavefunction was then relaxed over the two fine elements 
separately. In all, this process involved solving 3 linear 
systems of size n, so that the computational effort was 
1.5 times that of relaxing over each volume separately, 
rather than 4 times as in the case of simultaneous re- 
laxation. Despite this decrease in computational effort, 
rates of convergence were comparable to those obtained 
by relaxing over both volumes at once. The decrease 
in computational effort is more pronounced in higher di- 
mensions: using intergrid transfer to relax over a hyper- 
cube of 2^* elements increases costs by a factor of 1 -|- 2"'* 
relative to relaxing over each element separately, while re- 
laxing over the entire hypercube simultaneously increases 
costs by a factor of 2^*^. 



VI. CONCLUSIONS 

Efficiently and accurately propagating a time depen- 
dent wavefunction is a fundamental problem of compu- 



tational quantum mechanics, with applications ranging 
across many areas of physics. This paper has addressed 
this problem by developing an iterative procedure for re- 
laxing a trial wavefunction toward a variationally opti- 
mum, action minimizing solution. This relaxation proce- 
dure is trivially parallelizeable for problems involving a 
local Hamiltonian, and does not rely on the Hamiltonian 
being positive definite, making it well suited for incorpo- 
ration into multigrid relaxation methods. A local Fourier 
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FIG. 1. (Color online) Log ratio of the residual norm be- 
fore and after relaxation vs wavenumber (horizontal) and fre- 
quency (vertical) for trial wavefunction f{x,t) = g*{'=^+'^*)^ 
calculated using a single volume element. Columns show con- 
vergence for a given number of spatial and temporal basis 
functions, while rows show convergence for different values of 
the parameters k and i/ in Eq. 1301 a) k = 0, = 0, b) 
K = 0.1, = 0, c) K = 1.0, = 0, d) K = 10.0, u ^ 0, e) 
K = 0, u ^ 1.0, f) K = 0, u = -1.0, g) K = 1.0, u = 1.0, h) 
K — 1.0, I' — —1.0. Red areas correspond to ratios greater 
than 1.0, where the relaxation procedure does not converge. 

analysis shows that this procedure yields rapid rates of 
convergence over a wide range of parameters in the limit 
that the basis set defined over a particular volume ele- 
ment is sufficiently large. 

Although the discussion in this paper has focused on 
the specific problem of the time dependent Scrodinger 
equation, the analysis depends primarily on the existence 
of a local Hamiltonian. As this is a very common con- 
dition in physical problems, it may in the future be de- 
sirable to extend this analysis to other computationally 
difficult problems. 
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FIG. 2. (Color online) Log ratio of the magnitude of the ac- 
cumulated action before and after relaxation vs wavenumber 
(horizontal) and frequency (vertical) for trial wavefunction 
f{x,t) — e*(*^+'^*\ calculated using a single volume element. 
Columns show convergence for a given number of spatial and 
temporal basis functions, while rows show convergence for dif- 
ferent values of the parameters k and f in Eq. 1301 a) k = 0, 
ly = 0, h) K = 0.1, u = 0, c) K ^ 1.0, u ^ 0, d) K ^ 10.0, 
u = 0, e) K ^ 0, = 1.0, f) K = 0, = -1.0, g) K = 1.0, 
v = 1.0, h) K = 1.0, v — —1.0. Red areas correspond to 
ratios greater than 1.0, where the relaxation procedure does 
not converge. 



9 




FIG. 3. (Color online) Log ratio of the residual norm be- 
fore and after relaxation vs wavenumber (horizontal) and fre- 
quency (vertical) for trial wavefunction f{x,t) = g»{'=^+'^*)^ 
calculated using two adjacent volume elements. Columns 
show convergence for a given number of spatial and tempo- 
ral basis functions, while rows show convergence for different 
values of the parameters k and f in Eq. 1301 a) k = 0, = 0, 
b) K = 0.1, = 0, c) K = 1.0, u = 0, d) K = 10.0, u = 0, 
e) K^O, u ^ 1.0, f) K = 0, ;/ = -1.0, g) k = 1.0, u = 1.0, 
h) Mi — 1.0, v — —1.0. Red areas correspond to ratios greater 
than 1.0, where the relaxation procedure does not converge. 
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FIG. 4. (Color online) Log ratio of the magnitude of the ac- 
cumulated action before and after relaxation vs wavenumber 
(horizontal) and frequency (vertical) for trial wavefunction 
f{x,t) — e'(*'^+'^*\ calculated using a two adjacent volume 
elements. Columns show convergence for a given number of 
spatial and temporal basis functions, while rows show conver- 
gence for different values of the parameters k and f in Eq. 
|30] a) K = 0, 1/ = 0, b) K = 0.1, = 0, c) K = 1.0, = 0, d) 
K = 10.0, u = 0, e) K = 0, u = 1.0, f) K = 0, z/ = -1.0, g) 
K = 1.0, v = 1.0, h) hi — 1.0, v = —1.0. Red areas correspond 
to ratios greater than 1.0, where the relaxation procedure does 
not converge. 
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FIG. 5. (Color online) Log ratio of the residual norm be- 
fore and after relaxation vs wavenumber (horizontal) and fre- 
quency (vertical) for trial wavefunction f{x,i) — ^ 
calculated for two adjacent elements using the intergrid trans- 
fer method. Columns show convergence for a given number 
of spatial and temporal basis functions, while rows show con- 
vergence for different values of the parameters k and f in Eq. 
|30] a) K = 0, 1/ = 0, b) K = 0.1, = 0, c) K = 1.0, = 0, d) 
K = 10.0, u = 0, e) K ^ 0, u = 1.0, f) K = 0, u = -1.0, g) 
K — 1.0, v — 1.0, h) K = 1.0, v — —1.0. Red areas correspond 
to ratios greater than 1.0, where the relaxation procedure does 
not converge. 




FIG. 6. (Color online) Log ratio of the magnitude of the ac- 
cumulated action before and after relaxation vs wavenumber 
(horizontal) and frequency (vertical) for trial wavefunction 
f{x,t) = g»('=^+"')^ calculated for two adjacent elements us- 
ing the intergrid transfer method. Columns show convergence 
for a given number of spatial and temporal basis functions, 
while rows show convergence for different values of the pa- 
rameters K and u in Eq. 1301 a.) n — 0, i/ = 0, h) k = 0.1, 

V ^ 0, c) K = 1.0, 1/ = 0, d) K = 10.0, 1/ = 0, e) K = 0, 
ly = 1.0, f) K = 0, i/ = -1.0, g) «: = 1.0, u = 1.0, h) «: = 1.0, 

V = —1.0. Red areas correspond to ratios greater than 1.0, 
where the relaxation procedure does not converge. 



